

#include "matrix.h"

#include <stdio.h>

void test_basic_matrix();
void test_matrix_transform();
void test_detninv();
void test_addnsub();
void test_row_swap();
void test_gaussian_elim();

vm_matrix *m1;
vm_matrix *m2;
vm_matrix *system;

int main(int argc, char **argv){

  // First load the matrices to do testing on.
  fprintf(stderr, "Loading matrices...\n");
  m1 = vm_read_matrix("../data/test_matrix1.mtx");
  m2 = vm_read_matrix("../data/test_matrix2.mtx");
  system = vm_read_matrix("../data/test_system.mtx");
  fprintf(stderr, "Matrices loaded.\n");

  vm_print_matrix(m1, " ");
  vm_print_matrix(m2, " ");

  /* Now call test_basic_matrix() function. */
  //test_basic_matrix();
  /* Now test the transformations... */
  //test_matrix_transform();
  /* Do some tests with the determinant. */
  test_detninv();
  //test_addnsub();
  /* Test row swapping operations. */
  //test_row_swap();
  /* Test guassian elim */
  //test_gaussian_elim();
  
  /*vm_print_matrix(vm_matrix_mult(vm_read_matrix("../data/blah.mtx"), 
				 vm_read_matrix("../data/det2.mtx")));

    vm_print_matrix(vm_matrix_mult(m2, m1));*/


}


void test_basic_matrix(){

  vm_matrix_rc rc_row;
  vm_matrix_rc rc_col;

  vm_get_matrix_row(m1, 1, &rc_row);
  vm_get_matrix_col(m2, 1, &rc_col);
  vm_print_matrix_rc(&rc_row);
  vm_print_matrix_rc(&rc_col);

  double dot;
  if (!vm_matrix_dot(&rc_col, &rc_row, &dot)){
    printf("Error occured taking dot product.\n");
  } else {
    printf("Dot product: %f\n", dot);
  }

  vm_matrix *mult = vm_matrix_mult(m1, m2);
  if (mult !=  NULL){
    vm_print_matrix(mult, " ");
  } else {
    printf("Error occured during multiplication of matrices.\n");
  }

  vm_matrix *mtx = vm_read_matrix("../data/transpose.mtx");
  vm_matrix *mtx_t = vm_transpose(mtx);

  vm_print_matrix(mtx, " ");
  vm_print_matrix(mtx_t, " ");

  printf("Minor matrices:\n");
  int i, j;
  for ( i = 0; i < mtx->rows; i++)
    for ( j = 0; j < mtx->columns; j++){
      //vm_matrix *m = vm_minor_matrix(mtx, i, j);
      //vm_print_matrix(m);
      double cof = vm_cofactor(mtx, i, j);
      printf("Cofactor: %f\n", cof);
    }


}

void test_matrix_transform(){

  /* First load the transformation matrix. */
  vm_matrix *trans_m = vm_read_matrix("../data/transform3.mtx");
  if (trans_m == NULL) 
    fprintf(stderr, "Transformation matrix is null. Whoops\n");
  
  printf("Transformation vector:\n");
  vm_print_matrix(trans_m, " ");

  FILE *vectors = fopen("../data/transform.vec", "r");

  /* Now load and trasform vectors until the transform vectors run out. */
  while (1){
    vm_vector cur_vec;
    if (!vm_read(vectors, &cur_vec)){
      printf("Done reading vectors.\n");
      break;
    } else {
      vm_vector transformed;
      vm_print(cur_vec);
      printf(" ==> ");
      vm_transform(&cur_vec, trans_m, &transformed);
      vm_println(transformed);
    }
  }

  printf("Determinant of transformation matrix: ");

  double det;
  if (!vm_determinant(trans_m, &det)){
    printf("Error.\n");
  }
  printf("%f\n", det);

}

void test_detninv(){

  vm_matrix *det = vm_read_matrix("../data/det3.mtx");

  printf("----------------------------\nMatrix:\n");
  vm_print_matrix(det, " ");

  /*
  printf("Adjoint of the matrix:\n");
  vm_matrix *adj = vm_adjoint(det);
  vm_print_matrix(adj, " ");
  */

  //vm_rref(det);
  //vm_print_matrix(rref, " ");

  printf("Inverse of the matrix:\n");
  vm_matrix *inv = vm_inverse(det);
  vm_print_matrix(inv, " ");

  vm_matrix *test = vm_matrix_mult(inv, det);
  printf("A * inv(A):\n");
  vm_print_matrix(test, " ");


  /*
  printf("Determinant of the matrix:\n");
  double deter;
  int ok = vm_determinant(det, &deter);
  printf("Determinant: %f\n", deter);
  deter = _vm_rec_nbyn_det(det);
  printf("Determinant: %f\n", deter);
  */

}

void test_addnsub(){

  printf("----------------------------\n");
  printf("Adding two matrices:\n");
  vm_matrix *res = vm_matrix_add(m1, m2);
  if (res == NULL){
    printf("Error detected.\n");
  } else {
    vm_print_matrix(res, " ");
  }
}

void test_row_swap(){

  printf("Testing the matrix row swapping functionality.\n");

  if (vm_swap_rows(m1, 0, 2)){
    vm_print_matrix(m1, " ");
  } else {
    printf("Could not swap rows 0 and 2 in m1 for some reason.\n");
  }

  if (vm_swap_rows(m2, 1, 2)){
    vm_print_matrix(m2, " ");
  } else {
    printf("Could not swap rows 1 and 2 in m2 for some reason.\n");
  }


}

void test_gaussian_elim(){

  printf("Testing guassian elimination routines:\n");
  printf("  Starting matrix:\n");
  vm_print_matrix(system, " ");
  
  vm_row_echelon_form(system);
  printf("New row-echelon form matrix:\n");
  vm_print_matrix(system, " ");

  printf("Reducing the matrix to reduced row echelon form:\n");
  vm_rref(system);
  vm_print_matrix(system, " ");

}
